Heat transport and phonon localization in mass-disordered harmonic crystals 
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We investigate the steady state heat current in two and three dimensional disordered harmonic 
crystals in a slab geometry, connected at the boundaries to stochastic white noise heat baths at 
different temperatures. The disorder causes short wavelength phonon modes to be localized so the 
heat current in this system is carried by the extended phonon modes which can be either diffusive 
or ballistic. Using ideas both from localization theory and from kinetic theory we estimate the 
contribution of various modes to the heat current and from this we obtain the asymptotic system 
size dependence of the current. These estimates are compared with results obtained from a numerical 
evaluation of an exact formula for the current, given in terms of a frequency transmission function, 
as well as from direct nonequilibrium simulations. These yield a strong dependence of the heat ffux 
on boundary conditions. Our analytical arguments show that for realistic boundary conditions the 
conductivity is finite in three dimensions but we are not able to verify this numerically, except in 
the case where the system is subjected to an external pinning potential. This case is closely related 
to the problem of localization of electrons in a random potential and here we numerically verify that 
the pinned three dimensional system satisfies Fourier's law while the two dimensional system is a 
heat insulator. We also investigate the inverse participation ratio of different normal modes. 

PACS numbers: 



I. INTRODUCTION 

Energy transport in dielectric crystals at low temperatures 
is limited by isotope mass disorder and by anharmonicities. 
For a crystal coupled to thermal heat baths, in a slab geom- 
etry, the total energy transported may in addition depend 
on how the bath is coupled to the boundaries of the bulk. In 
this paper we will ignore anharmonicities but study the two 
other effects in considerable depth. 

Based on empirical evidence, one generically expects the 
validity of Fourier's law, i.e. , for a slab in contact with heat 
reservoirs at different temperatures the average energy cur- 
rent, J, should be proportional to 1/N, with N being the slab 
length. There have been many attempts to derive Fourier's 
law from microscopic dynamics. Using very reasonable phys- 
ical assumptions of local equilibrium and the notion of the 
mean free path traveled by phonons between collisions yields 
a heuristic derivation of Fourier's law. There is however no 
fully convincing derivation. There have also been many com- 
puter simulations of the heat flux in the nonequilibrium sta- 
tionary states of anharmonic crystals kept in contact with 
thermal reservoirs at different temperatures, and theoretical 
analysis based on the Green-Kubo formalism [11 [H E] ■ These 
studies suggest (but there is no proof despite some claims) 
that Fourier' law is not valid for one and two dimensional 
systems, even in the presence of anharmonic interactions, 
unless the system is also subjected to an external substrate 
pinning potential. Generically it is found that, for a system 
in contact with heat reservoirs at different fixed tempera- 
tures, the heat current density J scales anomalously with 
system length N as 



with /i 7^ 1 . The effective thermal conductivity behaves then 
as K ^ where a = I — fi. For two dimensional systems 
there are some analytic studies which suggest k ^ ln{N). 



Recent experiments on heat conduction in nanotubes and 
graphene flakes have reported observations which indicate 
such divergence of k with system size [H |S] . 

For heat conduction in the ordered harmonic crystal there 
are exact results from which one has /i = in all dimensions 
[6] 13 • Heat conduction in a disordered harmonic crystal 
will be affected by Anderson localization [8] and by phonon 
scattering. 

In this paper we report results of heat conduction studies 
in 2D and 3D disordered harmonic lattices with scalar dis- 
placements, connected to heat baths modeled by Langevin 
equations with white noise. We pay particular attention to 
the interplay between localization effects, boundary effects, 
and the role of long wavelength modes. The steady state 
heat current is given exactly as an integral over all frequen- 
cies of a phonon transmission coefficient. Using this for- 
mula and heuristic arguments, based on localization theory 
and kinetic theory results, we estimate the system size de- 
pendence of the current. The main idea behind our argu- 
ments is that the phonon states can be classified as ballistic 
modes, diffusive modes and localized modes. The classifi- 
cation refers both to the character of the eigenfunctions as 
well as to their transmission properties. Ballistic modes are 
spatially extended and approximately periodic; their trans- 
mission is independent of system size. Diffusive modes are 
extended but non-periodic and their transmission decays as 
1/7V. For localized modes transmission decays exponentially 
with N. In the context of kinetic theory calculations, the bal- 
listic modes are the low frequency modes with phonon mean 

free path ixi^) ^ N, and their contribution to the current 
leads to divergence of the thermal conductivity. Here we will 
carefully examine the effect of boundary conditions on these 
modes. 

Numerically we use two different approaches to study the 
nonequilibrium stationary state. The first is a numerical one 
which relies on the result that the current can be expressed 
in terms of a transmission coefficient. This transmission co- 
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efficient can be written in terms of phonon Green's functions 
and we implement efficient numerical schemes to evaluate 
this. The second approach is through direct nonequilibrium 
simulations of the Langevin equations of motion and finding 
the steady state current and temperature profiles. We have 
also studied properties of the isolated system, i.e., of the 
disordered lattice without coupling to heat baths and looked 
at the normal mode frequency spectrum and the wavefunc- 
tions. One measure of the degree of localization of the normal 
modes of the isolated system is the so-called inverse partici- 
pation ratio [IPR, defined in Eq. ( [To| below]. We have car- 
ried out studies of the IPR and linked these with the results 
from the transmission study. 

Phonon localization: This is closely related to the elec- 
tron localization problem. The effect of localization on linear 
waves in disordered media has been most extensively studied 
in the context of the Schrodinger equation for non-interacting 
electrons moving in a disordered potential. Looking at the 
eigenstates and eigenfunctions of the isolated system of a 
single electron in a disordered potential one finds that, in 
contrast to the spatially extended Bloch states in periodic 
potentials, there are now many eigenfunctions which are ex- 
ponentially localized in space. It was argued by Mott and 
Twose [S] and by Borland [10', and proven rigorously by 
Goldsheid et al. jjllj, that in one dimension {ID) all states 
are exponentially localized. In two dimensions (2D) there is 
no proof but it is believed that again all states are localized. 
In three dimensions (SD) there is expected to be a transition 
from extended to localized states as the energy is moved to- 
wards the band edges [T^. The transition from extended to 
localized states, which occurs when the disorder is increased, 
changes the system from a conductor to an insulator. The 
connection between localization and heat transport in a crys- 
tal is complicated by the fact that phonons of all frequencies 
can contribute to energy transmission across the system. In 
particular account has to be taken of the fact that low fre- 
quency phonon modes are only weakly affected by disorder 
and always remain extended. The heat current carried by 
a mode which is localized on a length scale £, decays with 
system length TV as e~^l^ . This I depends on the phonon 
frequency and low frequency modes for which I '-^ N will 
therefore be carriers of the heat current. The net current 
then depends on the nature of these low frequency modes 
and their scattering due to boundary conditions (BCs). 

A renormalization group study of phonon localization in 
a continuum vector displacement model was carried out by 
John et al. |13j . They found that much of the predic- 
tions of the scaling theory of localization for electrons carry 
over to the phonon case. Specifically they showed that in 
one and two dimensions all non-zero frequency phonons are 
localized with the low frequency localization length diverg- 
ing as i ^ uj^^ and ^ e'^^^ , respectively (where c > 
is some constant). This means that in 11? all modes with 

Lj ^ uj^ = N^^/"^ are localized while in 2D all modes with 

Lo ~ bj^ = [log(A^)]~^/^ are localized. In iD the prediction is 
that there is an uj^ independent of TV above which all modes 
are localized. However this study does not make any state- 
ments on the system size dependence of the conductivity. 

Kinetic theory: If one considers the low frequency ex- 
tended phonons, then the effect of disorder is weak and in 
dimensions d > \ one expects that localization effects can 



be neglected and kinetic theory should be able to provide an 
accurate description. In this case one can think of Rayleigh 
scattering of phonons. This gives an effective mean free path 
~ u}^^'^^^^ [see appendix [a] , for dimensions d > 1, 
and a diffusion constant D{uj) = viKito) where v, the sound 
velocity, can be taken to be a constant. For a finite system of 
linear dimension N we have D{uj) = vN for uj ^ ^ 
Kinetic theory then predicts 

K= dup{uj)D{u) , (2) 

where p(uj) ^ ui'^^^ is the density of states and we get 
K ^ implying fj, = d/{d+l). The divergence of the 

phonon mean free path at low frequencies and the resulting 
divergence of the thermal conductivity of a disordered har- 
monic crystal has been discussed in the literature and it has 
been argued that anharmonicity is necessary to make k finite 

Simulation results: There have been only few simula- 
tion studies of heat conduction in three dimensional disor- 
dered systems and none have been definitive concerning the 
validity of Fourier's law |16L I17j . In two dimension a diverg- 
ing thermal conductivity was reported in [18] . Some other 
studies have also looked at heat conduction in glassy systems 
at low temperatures where the harmonic approximation was 
used [ini [2U] but these did not address the questions of A^- 
dependence of k and the validity of Fourier's law. 

In ID it is well known from rigorous results and numer- 
ical studies that a ^ and its value is strongly dependent 
on boundary conditions [HI [21 [13 [Ml [23 For fixed 
BCs one has a = —1/2 while for free BCs, a = 1/2. The 
precise definitions of the different BCs will be given later. 
Here we explain the different physical situations they cor- 
respond to. If we model the heat reservoirs themselves by 
infinite ordered harmonic crystals then Langevin type equa- 
tions for the system 12 7 j are obtained on eliminating the bath 
degrees of freedom. The two different BCs then emerge nat- 
urally. Fixed BCs correspond to reservoirs with properties 
different from the system (e.g. different spring constants) 
and in this case one finds that effectively the particles at the 
boundaries (those coupled to reservoirs) experience an addi- 
tional harmonic pinning potential. Free BCs correspond to 
the case where the reservoir is simply an extension of the sys- 
tem (without disorder) and in this case the end particles are 
unpinned. Free BCs have been studied in the literature in 
the context of heat conduction in one dimensional chains |22j 
and in studies on nanotubes [28J . In this paper we study 
lattices with both fixed and free BCs although we think that 
fixed BCs are more realistic. 

In the presence of an external pinning potential low fre- 
quency modes are suppressed, hence one expects qualitative 
differences in transport properties. The pinned system has 
often been used as a model system to study the validity of 
Fourier's law. It has no translational invaraince and is thus 
more closely related to the problem of electrons moving in a 
random potential. Here we consider systems with and with- 
out external pinning potentials. 

The rest of the paper is organized as follows. In Sec. (|ll]) 
we define the specific model studied by us and present some 
general results for heat conduction in harmonic Hamiltonian 
systems connected to Langevin baths. We also give some 
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details of the numerical and simulation methods used in the 
paper. The transfer matrix approach used in evaluating the 
phonon transmission function is explained in AppendixIS] In 



Sec. ( III I a brief review of results for the one dimensional case 



and the heuristic arguments for the higher dimensional cases 
are given. In Sec. (IV I we present results from both the nu- 



merical approach and from nonequilibrium simulations. The 
main results presented are for transmission functions, IPRs of 
normal modes and the system size dependence of the current 
in two and three dimensional disordered harmonic lattices. 
Along the way we also present results for the density of states 
p{uj). Finally we conclude with a discussion in Sec. ([v|. 



II. MODELS AND METHODS 

For simplicity we consider only the case where longitudi- 
nal and transverse vibration modes are decoupled and hence 
we can describe the displacement at each site by a scalar 
variable. Also we restrict our study to d-dimensional hyper- 
cubic lattices. Let us denote the lattice points by the vector 
n = {ni, 712, n^} with = 1,2, ...,iV. The displacement 
of a particle at the lattice site n is given by Xn- In the har- 
monic approximation the system Hamiltonian is given by 



H 



^ 2 

n,e 



where e refers to the 2d nearest neighbors of any site and 
we impose boundary conditions which will be specified later. 
We have also included an external pinning harmonic poten- 
tial with spring constant fco, which we will sometimes set 
equal to zero. We consider a binary mass disordered crystal. 
Specifically we set the masses of exactly half the particles at 
randomly chosen sites to be ?7i — A and the rest to be m -I- A. 
Thus A gives a measure of the disorder. 

We couple all the particles at ni = 1 and nx — Nio heat 
reservoirs, at temperatures and respectively, and use 
periodic boundary conditions in the other (d — 1) directions. 
The heat conduction takes place along the 1 direction. Each 
layer with constant ni consists of N' = iV'^"-'^ particles. The 
heat baths are modeled by white noise Langevin equations 
of motion for the particles coupled to the baths. Using the 
notation n — (ni, n'), the equations of motion are given by: 

^nXn = y ' k[Xn 2;n-)-e) koXn + (5ni,l( T-^n 
e 

-I- 77^, - k'Xn) + SmM{~lXn + Vn' - k'Xn) , (4) 

where the dissipative and noise terms are related by the usual 
fluctuation dissipation relations 



(5) 



The particles at the surfaces ni = I, N experience additional 
harmonic pinning potentials with spring constants k' arising 
from coupling to the heat reservoirs. We consider two kinds 
of boundary conditions at the surfaces connected to reser- 
voirs: (i) fixed BCs fc' > and (ii) free BCs k' = 0. As 
discussed in the introduction fixed BCs correspond to reser- 
voirs with properties different from the system while free 



reservoirs correspond to the case where the reservoir is re- 
ally an extension of the system but without disorder. For the 
pinned case we only consider fixed BC. A schematic of the 
models and the different boundary conditions that we have 
studied is given in Fig. ([T]). 

Henceforth we will use dimensionless variables: force- 
constants are measured in units of k, masses in units of 
the average mass m, time in units of the inverse frequency 
fi^^ = (m/fc)^/^, displacements are in units of the lattice 
spacing a, friction constant 7 is in units of mil, and finally 
temperature is measured in units of fha^fl^ /ks- 

Driven by the reservoirs at two different temperatures 
and Tji the system reaches a nonequilibrium steady state. 
Our main interest will be in the steady state heat current in 
the system. Given the Langevin equations of motion Eq Q, 
one can find a formal general expression for the current. Let 
us denote by X a column vector with TV^ elements consist- 
ing of the displacements at all lattice sites. Similarly let 
X represent the vector for velocities at all sites. Then we 
can write the Hamiltonian in Eq. ([3| in the compact form 
H = ^X^AiX + ^X'^VX , which defines the diagonal mass 
matrix M and the force constant matrix V. With this nota- 
tion we have the following form for the steady state current 
per bond from the left to the right reservoir [531 127] '■ 



J = 



AT 

4:7:N' 



(iw7jv(w) , 



where 



(6) 



(7) 



and AT = Tl — Tjf . The , 5^ represent terms arising 
from the coupling to the left and right baths respectively, 
and Tl ji — Im[S'^ j^]. The specific form of for our 

system described by Eqs. Q is given in Appendix |b] The 
matrix g~^(uj) can be identified as the phonon Green's func- 
tion of the system with self-energy corrections due to the 
baths 127]. The integrand in Eq. ^ Tn{^) can be thought 
of as the transmission coefficient of phonons at frequency 
Lo from the left to the right reservoir. It will vanish, when 
A'^ — > 00 , at values of lo for which the disorder averaged den- 
sity of states is zero. Note that due to the harmonic nature 
of the forces the dependence of the heat flux on the reservoir 
temperatures enters only through the term AT in Eq. (|6|. 
The above expression for the current is of the Landauer form 
and has been derived using various other approaches such as 
scattering theory |30l I31j and the nonequilibrium Green's 
function formalism [32 [33] . 

Numerical approach: In Appendix [B] we describe how 
7jv can be expressed in a form amenable to accurate numer- 
ical computation. The system sizes we study are sufficiently 
large so that Tn{uj) has appreciable values only within the 
range of frequencies of normal modes of the isolated system, 
i.e., corresponding to 7 = in Eq. Q. Outside this range we 
flnd that the transmission rapidly goes to zero. By perform- 
ing a discrete sum over the transmitting range of frequencies 
we do the integration in Eq. ([6]) to obtain the heat current 
density J7. In evaluating the discrete sum over lo, step sizes 
of Suj = 0.01 — 0.0001 are used and we verified convergence in 
most cases. With our choice of units we have k = l,fh = 1 
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and we fixed AT = 1. Different values of the mass vari- 
ance A and the on-site spring constant kg were studied for 
two and three dimensional lattices of different sizes. It is ex- 
pected that the value of the exponent fj, will not depend on 
7 and in our calculations we mostly set 7=1, except when 
otherwise specified. 

Simulation approach: The simulations of Eq. Q are 
performed using a velocity- Ver let scheme as given in |34j . 
The current and temperature profiles in the system are ob- 
tained from the following time averages in the nonequilib- 
rium steady state: 



Ji 

Jn+1 

Tn 



Tl - ™(l,n')(i(l,n')/ 



n' \ I 



-y 



7 



(JV,n') 



n = 2,3,...,iV , 

Tr - '71(Ar,n')(i(JV,n')) 



We then obtained the average current J = 
^^n=i + !)• I'^ the steady state one has 

Jn = J for all n and stationarity can be tested by checking 
how accurately this is satisfied. We chose a step size of 
Ai = 0.005 and equilibrated the system for over 10* time 
steps. Current and temperature profiles were obtained by 
averaging over another 10* time steps. The parameters 
Tl = 2.0, Tr — 1.0 are kept fixed and different values of 
the mass variance A and the on-site spring constant kg are 
simulated. 

The value of Tn{lu), J and r„ depend, of course, on the 
particular disorder realization. Mostly we will here be inter- 
ested in disorder averages of these quantities which we will 
denote by [T], J = \J\ and [T„]. We also define the disorder 
averaged transmission per bond with the notation 

TH = ^[rH] . 

Numerical analysis of eigenmodes and eigenfunc- 
tions: We have studied the properties of the normal modes 
of the disordered harmonic lattices in the absence of cou- 
pling to reservoirs, again with both free and fixed boundary 
conditions. The d-dimensional lattice has p = 1,2, ...jA^'' 
normal modes and we denote the displacement field corre- 
sponding to the p**^ mode by a„(p) and the corresponding 
eigenvalue by Wp. The normal mode equation corresponding 
to the Hamiltonian in Eq. (|3| is given by: 



^nWpan = (2d -f fco)an - ^ 1 



n+e 



(8) 



where the an satisfy appropriate boundary conditions. In- 

1/2 

troducing variables V'n(p) = in(p), "^-a = (2(i -I- ko)/iin-n 
and tn,i = l/(mnmi)^/^ for nearest neighbour sites n, 1 the 
above equation transforms to the following form: 



(9) 



This has the usual structure of an eigenvalue equation for a 
single electron moving in a d-dimensional lattice correspond- 
ing to a tight-binding Hamiltonian with nearest neighbour 
hopping tn,i and on-site energies w„. Note that t^A and 
are correlated random variables, hence the disorder-energy 
diagram might differ considerably from a single band Ander- 
son tight-binding model. 

We have numerically evaluated all eigenvalues and eigen- 
states of the above equation for finite cubic lattices of size 
upto = 64 in 2D and = 16 in 3D. One measure of the 
degree of localization of a given mode is the inverse partici- 
pation ratio (IPR) defined as follows: 



P-1 = 



En 4 



(10) 



For a completely localized state, i.e. o„ = (5n,noi 



p-i takes 

the value 1. On the other hand for a completely delocalized 
state, for which On = 7V~'^/2e*"-i where q is a wave vector, 
takes the value N~'^. We will present numerical results 
for the IPR calculated for all eigenstates of given disorder 
realizations, in both 2D and 3D. Finally we will show some 
results for the density of states, p{uj), of the disordered sys- 
tem defined by: 



CO) 



(11) 



The density of states of disordered binary mass harmonic 
crystals was studied numerically by Payton and Visscher in 
1967 [35] and reviewed by Dean in 1972 [36^. 



III. HEAT CONDUCTION IN DISORDERED 
HARMONIC CRYSTALS: GENERAL 
CONSIDERATIONS 

Let us first briefiy consider heat conduction in the one di- 
mensional disordered harmonic chain. This has been exten- 
sively studied and is well understood [?Tl|^[^i24. 25, 2^. 
The matrix formulation explained in the last section leads to 
a clear analytic understanding of the main results. The cur- 
rent is given by the general expression Eq. mj . From Eq. ( B5 1 
the transmission is given by Tn{uj) — 47^c[?|G^(w)p w 



tiere 

G~l^{uj) is now just a complex number. The disorder averaged 
transmission is given by Tn(lu) — [T/v(w)]. There are three 
observations that enable one to determine the asymptotic 
system size dependence of the curren t. Th es e are : 

(i) p(i.^) = [G^]^i given by Eqs. ( |B17| , \B19\ is a com- 
plex number which can be expressed in terms of the prod- 
uct of N random 2x2 matrices. Using Furstenberg's the- 
orem it can be shown that for almost all disorder realiza- 
tion, the large N behaviour of P^^-^) for fixed w > is 
|p(i'^) I ~ gbNuj ^ ^iiQYe 5 > is a constant. This is to be un- 
derstood in the sense that limAr^oo(l/-^) log IP""^'^-* | ~ 

for to ^ 0. Since Tn{uj) - |p(i.w)|-2 ^ ^-2bN>^^ ^ ^j^-g 
plies that transmission is significant only for low frequencies 
u! ^ uJc{N) ^ XjN'^l'^. The current is therefore dominated 
by the small uj behaviour of Tj'^^lu). 

(ii) The second observation made in [23] is that the trans- 
mission for Lu < u!c{N) is ballistic in the sense that Tpf{uj) is 
insensitive to the disorder. 



5 



(iii) The final important observation is that the form of 
the prefactors of e~^^" in Tjv(a') for uj < u}c{N) depends 
strongly on boundary conditions and bath properties |25l 
25] , For the white noise Langevin baths one finds Tn{uj) ~ 
for fixed BC and Tn{oj) ~ w°e-^^'^' for free BC 
[25] , This difference arises because of the scattering of long 
wavelength modes by the boundary pinning potentials. 

In Fig. ([2]) we plot numerical results showing Tjf{uj) for 
the ID binary mass-disordered lattice with both fixed and 
free boundary conditions. One can clearly see the two fea- 
tures discussed above namely (i) dependence of frequency 
cut-off on system size and (ii) dependence of form of Tj\j (w) 
on boundary conditions. Using the three observations made 
above it is easy to arrive at the conclusion that J ~ N~^/'^ 
for fixed BC and J ^ N~^^^ for free BC. In the presence of 
a pinning potential the low-frequency modes are suppressed 
and one obtains a heat insulator with J ~ e~'^^ , with c a 
constant p4] (see also |37j and references there). 

Higher dimensions. Let us try to extend the analysis 
of the \D case to higher dimensions. For this we will use 
inputs from both kinetic theory and the theory of phonon 
localization. The main point of our arguments involves the 
assumption that normal modes can be classified as ballistic, 
diffusive or localized. Using localization theory we determine 
the frequency region where states are localized. The lowest 
frequency states with lo ^ Q will be ballistic and we use 
kinetic theory to determine the fraction of extended states 
which are ballistic. We assume that at sufficiently low fre- 
quencies the effective disorder is always weak (even when the 
mass variance A is large) and one can still use kinetic theory. 
Corresponding to the three observations made above for the 
ID case we now make the following arguments: 

(i) From localization theory one expects all fixed non-zero 
frequency states in a 2D disordered system to be localized 
when the size of the system goes to infinity. As discussed 
in Sec. Q localization theory gives us a frequency cut-off 
uj^ = (lniV)~^/^ in 2D above which states are localized. In 
3-D one obtains a finite frequency cut-off uj^ independent of 
system size above which states are localized. 

(ii) For the unpinned case with finite N there will exist low 
frequency states below w^?', in both 2D and 31?, which are 
extended states. These states are either diffusive or ballis- 
tic. Ballistic modes are insensitive to the disorder and their 
transmission coefficient are almost the same as for the or- 
dered case. To find the frequency cut-off below which states 
are ballistic we use kinetic theory results (see Appendix |A]) . 
For the low-frequency extended states we expect kinetic the- 
ory to be reliable and this gives us a mean free path for 
phonons tK ~ 0;"^'^+^^. This means that for low frequen- 
cies to ^ u!^ ~ ]\[-'^/(d.+i) have £k{(^) > N and phonons 
transmit ballistically. We now proceed to calculate the con- 
tribution of these ballistic modes to the total current. This 
can be obtained by looking at the small uj form of Tn{u!) for 
the ordered lattice. 

(iii) For the ordered lattice Tn{uj) is typically a highly 
oscillatory function with the oscillations increasing with sys- 
tem size. An effective transmission coefficient in the N ^ 00 
limit can be obtained by considering the integrated trans- 
mission. This asymptotic effective low-frequency form of 
Tpf{uj), for the ordered lattice can be calculated using meth- 



ods described in 



T{co) 
T{u) 



and is given by: 

- uj'^+^ , fixed BC 

- uj'^-'^ , free BC . 



(12) 



the result being valid for d = 1, 2, 3 P8]. 

Using the above arguments we then get the ballistic con- 
tribution to the total current density (for the unpinned case) 



d+1 



duj bj' 



dTT) ' fixed BC, 



1 



]^d/(d+ 



-r , free BC . (13) 



We can now make predictions for the asymptotic system size 
dependence of total current density in two and three dimen- 
sions. 

Two dimensions: From localization theory one expects 
that all finite frequency modes uj ^ lj^ = {h\N)~^/'^ are 
localized and their contribution to the total current falls 
exponentially with system size. Our kinetic theory argu- 
ments show that the low frequency extended states with 
oj^ ^ UJ ^ oj^ are diffusive (where uj^ = N^^/'^) while 
the remaining modes with lo ^ uj^ are ballistic. The diffu- 
sive contribution to total current will then scale as Jdiff ^ 
{\aN)~^/'^N^^. The ballistic contribution depends on BCs 
and is given by Eq. (13). This gives Jbaii ~ A^~*/^ for fixed 
BC and Jbaii iV~" for free BC. Hence, adding all the 
different contributions, we conclude that asymptotically: 



J 



1 

[In Ny/'^N ' 
1 



fixed BC, d = 2 , 



free BC, d = 2. 



(14) 



In the presence of an onsite pinning potential at all sites the 
low frequency modes get cut off and all the remaining states 
are localized, hence we expect: 



J. 



-bN 



pinned , d — 2 



(15) 



where b is some positive constant. 

Three dimensions: In this case localization theory tells us 

that modes with ui ^ uj^ are localized and uj^ is indepen- 
dent of N. From kinetic theory we find that the extended 



states with lu^^ 



are diffusive (with = N~^^^) 



and those with uj ~ lu^ are ballistic. The contribution to 
current from diffusive modes scales as Jdiff ~ iV~^. The bal- 
listic contribution (from states with uj ~ N^^^^) is obtained 



from Eq. (13) and gives Jbaii ~ N^^^^ for fixed BC and 
Jbaii ^ for free BC. Hence, adding all contributions, 

we conclude that asymptotically: 



iV3/4 



fixed BC , d = 3 , 
free BC , d = 3 



(16) 



In the presence of an onsite pinning potential at all sites the 
low frequency modes get cut off and, since in this case the 
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pinned , d = 3. 



(17) 



remaining states form bands of diffusive and localized states, 
hence we expect: 

1 

N 

Thus in 3D both the unpinned lattice with fixed bound- 
ary conditions and the pinned lattice are expected to show 
Fourier type of behaviour as far as the system size depen- 
dence of the current is considered. 

Note that for free BC, the prediction for the current contri- 
bution from the ballistic part Jbaii ^ iV^'*/^'^"'"^) is identical 
to that from kinetic theory discussed in Sec. This agree- 
ment can be traced to the small ui form of T(uj) ^ lu'^~^ for 
free BC [see Eq. (12)] which is identical to the form of the 
density of states p(lu) used in kinetic theory. The typical 
form of density of states for ordered and disordered lattices 
in different dimensions is shown in Fig. ([s]) and we can see 
that the low frequency form is similar in both cases and has 
the expected behaviour. However it seems reasonable 

to expect that, since the transport current phonons are in- 
jected at the boundaries, in kinetic theory one needs to use 
the local density of states evaluated at the boundaries. For 
fixed BC this will then give rise to an extra factor of 
(from the squared wavefunction) and then the kinetic theory 
prediction matches with those given above. 

We note that the density of states in Fig. ^ show appar- 
ent gaps in the middle ranges of lu for d = 2,3. These might 
be expected to disappear when the size of the system goes to 
infinity when there should be large regions containing only 
masses of one type [231 HI] • These regions will however be 
rare. In Fig. Q we show plots of the density of states for 
the ordered and disordered harmonic lattices in the presence 
of pinning. In this case the gaps in the spectrum are more 
pronounced and, for large enough values of kg and A, may 
be present even in the thermodynamic limit. 



IV. RESULTS FROM NUMERICS AND 
SIMULATIONS 

We now present the numerical and simulation results for 
transmission coefficients, heat current density, temperature 
profiles and IPRs for the disordered harmonic lattice in var- 
ious dimensions. The numerical scheme for calculating J is 
both faster and more accurate than nonequilibrium simula- 
tions. Especially, for strong disorder, equilibration times in 
nonequilibrium simulations become very large and in such 
cases only the numerical method can be used. However we 
also show some nonequilibrium simulation results. Their 
almost perfect agreement with the numerical results pro- 
vides additional confidence in the accuracy of our results. 
In Sec. (IV A I we give the results for the 2D lattice for the 



in the conducting direction = 1). One of the interesting 
questions here is as to how the three properties for the ID 



case discussed in Sec. (HI I get modified for the 2D case. 



unpinned case with both fixed and free boundary conditions 



and then for the pinned case. In Sec. (IVBl we present the 



results for the three dimensional case with and without sub- 
strate pinning potentials. 

A. Results in two dimensions 

In this section we consider N x N square lattices with pe- 
riodic BCs in the 1^ — 2 direction and either fixed or free BCs 



1. Disordered 2D lattice wittiout pinning 

Fixed BC : we have computed the transmission coefficients 
and the corresponding heat currents for different values of A 
and for system sizes from = 16 — 1024. The number of 
averages varied from over 100 samples for = 16 to about 
two samples for N — 1024. In Figs. ( 5|6|7 l we plot the 
disorder averaged transmission coefficient for three different 
disorder strengths, A = 0.95, A = 0.8 and A 0.2, for 
different system sizes. The corresponding plots of IPRs as 
a function of normal mode frequency ujp, for single disorder 
realizations, are also given. From the IPR plots we get an 
idea of the typical range of allowed normal mode frequencies 
and their degree of localization. Low IPR values which scale 
as 7V~^ imply extended states while large IPR values which 
do not change much with system size denote localized states. 
In Fig. ^ we also show typical plots of small IPR and large 
IPR wavefunctions. From Figs. ( 5|6|7 ) we make the following 
observations: 

(i) As expected we see significant transmission only over 
the range of frequencies with extended states. Thus in 
Fig. ([5]) for A = 0.95 we see that, while there are normal 
modes in the range w « (0 — 12), transmission is appreciable 
only in the range w ~ (0 — 1.5) and this is also roughly the 
range where the IPR data shows a iV~^ scaling behaviour. 
This can also be seen in Fig. ([6| where the inset shows the 
decay of T{uj) in the localized region. Unlike the ID case 
we see a very weak dependence on system size of the up- 
per frequency cut-off oj^ beyond which states are localized 
and transmission is negligible. As discussed earlier, local- 
ization theory predicts uj^ ^ (lniV)~^/^ but this may be 
difficult to observe numerically. The overall transmission 
function Tj^{uj) decreases with increasing system size, with 
T{u!) ~ 1/7V at higher freqencies and T{uj) ~ 7V° at the 
lowest frequencies. 

(ii) In Fig. ^ we have also plotted T{lu) for the ordered 
binary mass case and we note that over a range of small 
frequencies, T{uj) for the disordered case is very close to the 
curve for the ordered case, which means that these modes are 
ballistic. As expected from the arguments in Sec. (HI I we 
roughly find T{uj) u)^ at small frequencies. The remaining 
transmitting states are either diffusive (with a 1/A'^ scaling) 
or are in the cross-over regime between diffusive and ballistic 
and so do not have a simple scaling. 

We next look at the the integrated transmission which 
gives the net heat current. The system size dependence of 
the disorder averaged current J for different values of A is 
shown in Fig. ([8|. For the case A = 0.2, we also show simula- 
tion results and one can see that there is excellent agreement 
with the numerical results. For A = 0.2 we get an exponent 
/z « 0.6 which is close to the value obtained earlier in [TB] 
for a similar disorder strength. However with increasing dis- 
order we see that this value changes and seems to settle to 
around /i ~ 0.75. It seems reasonable to expect (though we 
have no rigorous arguments) that there is only one asymp- 
totic exponent and for small disorder one just needs to go 
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to very large system sizes to see the true value. In Fig. ([9| 
we show temperature profiles obtained from simulations for 
lattices of different sizes with A = 0.2. The jumps at the 
boundaries indicate that the asymptotic size limit has not 
yet been reached. This is consistent with our result that the 
exponent ^ obtained at A = 0.2 is different from what we 
believe is the correct asymptotic value (obtained at larger 
values of A). We do not have temperature plots at strong 
disorder where simulations are difficult. 



Thus contrary to the arguments in Sec. (Ill I which pre- 
dicted J - (ln7V)-i/2iv-i we find a much larger current 
scaling as J ^ Af~°-^^. It is possible that one needs to go to 
larger system sizes to see the correct scaling. 



Free BC: In this case from the arguments in Sec. (Ill I we 



expect ballistic states to contribute most significantly to the 
current density giving J ~ N~'^/^ . 

In Figs. ( 10|11 1 we plot the disorder averaged transmission 
coefficient for A = 0.8 and A = 0.2 for different system sizes. 
Qualitatively these results look very similar to those for fixed 
boundaries. However transmission is now significantly larger 
in the region of extended states. The behaviour at frequen- 
cies w ^ is also different and we now find T(a;) ~ a; in 
contrast to T{lli) • 



for fixed boundaries. From the plots 
of IPRs in Fig. (|10| we note that there is not much quali- 



tative difference with the fixed boundary plots except in the 
low frequency region (see below). 

The system size dependence of the disorder averaged cur- 
rent J for two different values of A is shown in Fig. (12 1. 
For A = 0.2 we get an exponent /i sa 0.5 while for the 
stronger disorder case A = 0.8 we see a different exponent 
/i « 0.6. Again we believe that the strong disorder value of 
/i — 0.6 is closer to the value of the true asymptotic expo- 
nent. This value is close to the expected /j, = 2/3 for free BC 
and significantly different from the value obtained for fixed 
BC {jjL « 0.75). Thus the dependence of the value of a on 
boundary conditions exists even in the 2D case. 

For the case of free BCs, we find that the values of T{uj) in 
the diffusive regime matches with those for fixed BCs but are 
completely different in the ballistic regime. This is seen in 
Fig. ( 13 1 where we plot the effective mean free path Zcff(w) — 
NT{uj)/w'^~^ in the low- frequency region [this is obtained 
by comparing Eq. |6| with the kinetic theory expression for 
conductivity Eq. (|2fJ. For free BC, ^cff is roughly consistent 
with the kinetic theory prediction l~g ^ + iji~^{uj) 

but the behaviour for fixed BC is very different. The inset of 



Fig. ( 13 1 plots Icff for the equal mass ordered case and we find 



that in the ballistic regime it is very close to the disordered 
case, an input that we used in the heuristic derivation. The 
numerical data also confirms that for small lu, T{lo) ~ uj for 
free BCs and as lo^ for fixed BCs. The transmission for fixed 
BC shows rapid oscillations which increase with system size, 
and arise from scattering and interference of waves at the 
interfaces. 



2. Disordered 2D lattice with pinning 

We now study the effect of introducing a harmonic pinning 
potential at all sites of the lattice. It is expected that this will 
cut off low frequency modes and hence one should see strong 
localization effects. The localization length ^ will decrease 



both with increasing A and increasing ko (i n \D heuristic 
arguments give (. ~ l/(A^A:o) ^Tj). In Figs. ( 14|15 1 we plot 
the transmission coefficients for two cases with on-site po- 
tentials ko — 10.0 and k^ — 2.0 respectively, and A — 0.4. 
We also plot the IPR in Fig. (14). Unlike in the unpinned 



case we now find that the transmission coefficients are much 
smaller and fall more rapidly with system size. 

From the plot of we find that for all the modes, the 
value of P^^ does not change much with system size which 
implies that all modes are localized. The allowed frequency 
bands correspond to the transmission bands. The two wave- 
functions plotted in Fig. ( 14 1 correspond to one relatively 
small and one large P"^ value and clearly show that both 
states are localized. 

The system size dependence of the integrated current is 
shown in Fig. (16) for the two parameter sets. The val- 
ues of /X « 1.6, 3.65 for the two sets indicate that at large 
enough length scales one will get a current falling exponen- 
tially with system size and hence we have an insulating phase. 
In Fig. (17) we plot the temperature profiles for the set with 
A = 0.4, ko = 10.0 . In this case it is difficult to obtain 
steady state temperature profiles from simulations for larger 
system sizes. The reason is that the temperature (unlike cur- 
rent) gets contributions from all modes (both localized and 
extended) and equilibrating the localized modes takes a long 
time. 



B. Results in three dimensions 

In this section we mostly consider N x N x N lattices 
with periodic boundary conditions in the = 2, 3 directions. 
Some results for N x N2 x N3 lattices with N2 = N3 < N 
will also be described. Preliminary results for the case of free 
BCs are given and indicate that there is no dependence of 
the exponent fj, on BCs. It is not clear to us whether this is 
related to the boundedness of the fiuctuations in Xn and the 
decay of the correlations between Xn and x\ (like |n — 1|^^) 
in d = 3 and their growth (with TV) in d < 3. 



1. Disordered 3D lattice without pinning 

Fixed BC : we have used both the numerical approach and 
simulations for sizes up to 32 x 32 x 32 for which we have data 
for T{u}). For larger systems the matrices become too big and 
we have not been able to use the numerical approach. Hence, 
for larger system sizes we have only performed simulations, 
including some on N x N2X N2 lattices. For these cases only 
the current J is obtained. The number of averages varies 
from over 100 samples for = 16 to two samples for N = 64. 
In Figs. ( 18|19 1 we plot the disorder averaged transmission 
coefficient for two different disorder strengths, A — 0.8 and 
A = 0.2, for different system sizes. The corresponding plots 
of IPRs as a function of normal mode frequency ujp , for single 
disorder realizations, are also given. From the IPR plots 
we get an idea of the typical range of allowed normal mode 
frequencies and their degree of localization. Low IPR values 
which scale as N^^ imply extended states while large IPR 
values which do not change much with system size denote 
localized states. 



8 



From Figs. ( 18|19 1 we make the following observations. 

(i) From the 31? data it is clear the effect of localization 
is weaker than in ID and 2D. Both for A ~ 0.2 and A = 
0.8 we find that there is transmission over almost the entire 
range of frequencies of the allowed normal modes. From the 
IPR plots we see that for A = 0.2 most states are extended 
except for a small region in the high frequency band-edge. 
For A = 0.8 the allowed modes form two bands and one 
finds significant transmission over almost the full range. At 
the band edges (except the one at uj — 0) there are again 
localized states. It also appears that there are some large IPR 
states interspersed within the high frequency band. As in the 
2D case and unlike the ID case, the frequency range over 
which transmission takes place does not change with system 
size, only the overall magnitude of transmission coefficient 
changes. 



(ii) The plot of NT{uj) in Fig. ( 18 ) shows the nature of the 



extended states. The high frequency band and a portion of 
the lower frequency band have the scaling T{u!) ~ and 
hence corresponds to diffusive states. In the lower-frequency 
band the fraction of diffusive states seems to be increasing 
with system size but it is difficult to verify the uj^ ~ TV"^/^ 
scaling. The ballistic nature of the low-frequency states is 



confirmed in Fig. ( 19 1 where we see that T^uj) for the binary- 



mass ordered and disordered lattices match for small uj [with 



a T{uj) 



dependence] 



In Fig. (20) we show the system size dependence of the dis- 



order averaged current density J for the two cases with weak 
disorder strength (A = 0.2) and strong disorder strength 
(A = 0.8). The results for cubic lattices of sizes up to = 32 
are from the numerical method while the results for larger 
sizes are from simulations. We find an exponent /i « 0.6 at 
small disorder and « 0.75 at large disorder strength. As in 
the 2D case here too we believe that at small disorder, the 
asymptotic system size limit will be reached at much larger 
system sizes and that the exponent obtained at large disor- 
der strength is probably close to the true asymptotic value. 
The value (/x = 0.75) does not agree with the prediction 
( J ~ N~^) made from the heuristic arguments in Sec. (Ill I. 
A study of larger system sizes is necessary to confirm whether 
or not the asymptotic size limit has been reached. 

The data point at TV = 128 for the set with A = 0.2 in 
Fig. ( 20 ) actually corresponds to a lattice of dimensions 128 x 



48 X 48 and we believe that the current value is very close 
to the expected fully 3D value. To see this point, we have 



plotted in Fig. (21 1 results from nonequilibrium simulations 



with N X N2X N2 lattices with N2 < N. 



Finally, in Fig. (22 1 we show temperature profiles (for sin- 



gle disorder realizations) obtained from simulations for lat- 
tices of different sizes and with A = 0.2. The jumps at the 
boundaries again indicate that the asymptotic system size 
limit has not been reached even at the largest size. 



Free BC: In this case from the arguments in Sec. (Ill I we 



expect ballistic states to contribute most significantly to the 
current density giving J ^ N~^^^. 



In Fig. (23) we plot the disorder averaged transmission co- 
efficient for A = 0.8 for different system sizes. The transmis- 
sion function is very close to that for the fixed boundary case 
except in the frequency region corresponding to non-diffusive 
states. At w — + we now expect, though it is hard to verify 
from the data, that T{uj) ^ uP' in contrast to T{u)) ~ uj'^ for 



fixed boundaries. 

The system size dependence of the disorder averaged cur- 
rent J for two different values of A is shown in Fig. ( 20 1 . We 



find that the current values are quite close to the fixed BC 
case and the exponent obtained at the largest system size 
studied for this case is /i « 0.71. This value is close to the 
expected ^ = 3/4 for free BC. 

We now compare the transmission coefficient for free 
and fixed BCs in the ballistic regime. This is plotted 
in Fig. (24) where we show the effective mean free path 
/eff(i^) = NTiui) /w'^~^ in the low-frequency region. As in the 
2D case we again find that for free BCs, leS is roughly consis- 
tent with the kinetic theory prediction 1~q ~ N~'^ 
and the behaviour for fixed BCs is very different. The inset 
of Fig. (24) plots Iqs for the equal mass ordered case and we 



find that in the ballistic regime it is very close to the dis- 
ordered case. The numerical data confirms the input in our 



theory on the form of T{u>) for small uj, 



i.e. Tiuj) ~ Lj for 
free BCs and as for fixed BCs. The transmission for fixed 
BC shows rapid oscillations which increase with system size, 
and arise from scattering and interference of waves at the 
interfaces. 



2. Disordered 3D lattice with pinning 

For the pinned case, we again use both the numerical 
method and simulations for sizes up to A^ = 32. For A^ = 64 
only nonequilibrium simulation results are reported. 

In Figs. ( 25|26 1 we plot the disorder averaged transmission 
coefficient for A = 0.2 and A = 0.8 with ko = 10.0. The 
corresponding IPRs P^^ and scaled IPRs N^P~^ are also 
shown. 

From the IPR plots we notice that the spectrum of the 3D 
disordered pinned chain has a similar interesting structure 
as in the 2D case with two bands and a gap which is seen at 
strong disorder. However unlike the 2D case where all states 
were localized, here the IPR data indicates that most states 
except those at the band edges are diffusive. We see localized 
states at the band edges and also there seem to be some lo- 
calized states interspersed among the extended states within 
the bands. The insets in Figs. ( 25|26 ) show that there is a 
reasonable N^^ scaling of the transmission data in most of 
the transmitting region. This is clearer at the larger system 
sizes. Thus, unlike the unpinned case where low frequency 
extended states were ballistic or super-diffusive, here we find 
that there is no transmittance at small (w — s- 0) frequencies 
and that all states are diffusive. 

From the above discussion we expect Fourier's law to be 
valid in the 31? pinned disordered lattice. The system size 
dependence of the disorder averaged current J for different 
disorder strengths is plotted in Fig. ( 27 1 . For all the param- 



eter sets the exponent obtained is close to fi — I correspond- 
ing to a finite conductivity and validity of Fourier's law. The 
temperature profiles plotted in Fig. ( 28 1 have small bound- 



ary temperature jumps and indicate that the asymptotic size 
limit has already been reached. 

One might expect that at very strong disorder, all states 
should become localized and then one should get a heat in- 
sulator. The parameter set corresponding to Fig. (26) corre- 



sponds to strong disorder and for this we still find a signif- 
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d = 2 


d = 3 


Analytical 


Numerical 


Analytical 


Numerical 


Pinned 


exp {-hN) 








Fixed 


iV-i(lnAr)-i/2 


iV-0.75 




jY-0.75 


Free 




^-0.6 


^-3/4 


^-0.71 



TABLE I: The table summarizes the main results of the paper. 
The numerical (and nonequilibrium simulation) results obtained 
in the paper are compared, in two and three dimensions, with 
the analytical predictions obtained from our heuristic arguments. 
The error bar for the numerically obtained exponent values is of 
the order ±0.02. This error is estimated from the errors in the 
last few points of the J-versus-A'^ data. NB: The system sizes used 
may well be far from asymptotic. 



icant fraction of extended states. Thus for the binary mass 
case it appears that there are always extended states. We 
have some results for the case with a continuous mass dis- 
tribution ( masses are chosen from a uniform distribution 
between 1 — A and 1 + A). In this case we find that the 
effect of disorder is stronger and the transmission at all fre- 
quencies is much reduced compared to the binary mass case. 
However we cannot see the exponential decrease in transmis- 
sion with system size and so it is not clear if an insulating 
behaviour is obtained. Further numerical studies are neces- 
sary to understand the asymptotic behaviour. 



V. DISCUSSION 

We have studied heat conduction in isotopically disor- 
dered harmonic lattices with scalar displacements in two and 
three dimensions. The main question addressed is the sys- 
tem size dependence of the heat current, which is computed 
using Green's function based numerical methods as well as 
nonequilibrium simulations. We have tried to understand the 
size dependence by looking at the phonon transmission func- 
tion Tioj) and examining the nature of the energy transport 
in different frequency regimes. We also described a heuris- 
tic analytical calculation based on localization theory and 
kinetic theory and compared their predictions with our nu- 
merical and simulation results. This comparison is summa- 
rized in Table Q. 

The most interesting findings of this work are: 
(i) For the unpinned system we find that in 2D there are 
a large number of localized modes for which phonon trans- 
mission is negligible. In ZD the number of localized modes 
is much smaller. The extended modes are either diffusive or 
ballistic. Our analytic arguments show that the contribution 
of ballistic modes to conduction is dependent on BCs and is 
strongly suppressed for the case of fixed BCs, the more re- 
alistic case. In ZD this leads to diffusive modes dominating 
for large system sizes and Fourier's law is satisfied. Thus 
a finite heat conductivity is obtained for the 3D disordered 
harmonic crystal without the need of invoking anharmonic- 
ity as is usually beheved to be necessary [TH [T^. This is 
similar to what one obtains when one adds stochasticity to 
the time evolution in the bulk as shown by [39' . Our numer- 
ical results verify the predictions for free BCs and we believe 
that much larger system sizes are necesary to verify the fixed 
BC results ( this is also the case in ID [25l [26]). 



(ii) In two dimensions the pinned disordered lattice shows 
clear evidence of localization and we obtain a heat insulator 
with exponential decay of current with system size. 

(iii) Our result for the 3D pinned disordered lattice provides 
the first microscopic verification of Fourier's law in a three 
dimensional system. For the binary mass distribution we do 
not see a transition to insulating behaviour with increasing 
disorder. For a continuous mass distribution we find that the 
current is much smaller (than the binary mass case with the 
same value of A) but it is not clear whether all states get 
localized and if an insulating phase exists. 
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APPENDIX A: KINETIC THEORY 

Kinetic theory becomes valid in the limit of small disorder. 
Its basic object is the Wigner function, /, which describes 
the phonon density in phase space and is governed by the 
transport equation 

^ /(r, k, t) + VLo{k) ■ Vrf{r, k, t) = C/(r, fc, t) . (Al) 



dt 

Here r E 



l'^ (boundary conditions could be imposed), k g 
[— tTjTt]'' is the wave number of the first BriouUin zone, uj is 
the dispersion relation of the constant mass harmonic crystal, 
and C is the collision operator. It acts only on wave numbers 
and is given by 



Cfik) = {2T:)-''+'Lo{kyA 



dk' 



5{i^{k)^uo{k')){f{k')-f{k)). (A2) 
We refer to |40] for a derivation. In the range of validity 



of (All, (A2| we can think of phonons as classical particles 
with energy lj and velocity \Juj{k). They are scattered by 
the impurities from k to dk' with the rate 

{2TT)-'^+^oj{kf ls?5{uj{k) - u{k'))dk' . (A3) 

Collisions are elastic. We distinguish 

(i) no pinning ■potential. Then for small k one has a'(fc) = 
|fc| and |Va;(fc)| — 1. From (A2) the total scattering rate 
behaves as This is the basis for the discussion in 
connection with Eq. ^ . 

(ii) pinning poten tial. In this case ijj{k) = wo + fc^ for small k. 

,2 



The prefactor in ( A2 1 can be replaced by Wg . The velocity is 
k and the scattering is isotropic with rate Thus the 

diffusion coefficient results as D(fc) = |fc|~'^+'* which vanishes 
as |fc| for d = 2, 3. Hence there is no contribution to the 
thermal conductivity from the small k modes. 



APPENDIX B: TRANSFER MATRIX APPROACH 

We now outline steps by which TJv can be expressed in 
forms which are amenable to accurate numerical evaluation. 
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We will give results whereby we express T/v in terms of prod- 
uct of random matrices. These are related to the Green's 
function and transfer matrix methods used earlier in the cal- 
culation of localization lengths in disordered electronic sys- 
tems [H]. Some related discussions for the phonon case can 
be found in |42j . For heat conduction in one dimensional dis- 
ordered chains, the transfer matrix approach has been shown 
to be very useful in obtaining analytic as well as accurate 
numerical results and here we study the extension of this to 
higher dimensions. 

The transmission coefficient is given by Tn(uj) = 
4Tr[lL{Lo)g+{uj)lR{oj)g-{Lj) where G+iui) = [-iv^M + V- 
—S^]~^ , Q~ — [G^]*, 1l,r = I'^iSt^Ri we now spec- 
ify the form of S'l ^ corresponding to the equations of motion 
in Eqs. Q. Note that we have transformed to dimensionless 
variables w ui/Vl,M M/ifi, V V/k, 7 —^ j/{rhfl) 
where — {k/fhy/^. We are considering heat conduction in 
the V — 1 direction of a d-dimensional lattice with particles 
on the layers ni — 1 and ni = N being connected to heat 
baths at temperatures and Tr respectively. The matri- 
ces S'^ and S'lf represent the coupling of the system to the 
left and right reservoirs respectively, and can be written as 
N X N block matrices where each block is a A^' x N' matrix. 
The block structures are as follows: 



We first introduce some notation. Let 3;('>'+"~i) with 
l<n<N — l + l denote a n x n tridiagonal block ma- 
trix whose diagonal entries are a/, aj+i, ...ai^n-i, where each 
ai is a, N' X N' matrix. The off-diagonal entries are given 
by — /. For an arbitrary block matrix A^'''™\ -^[i ™^ ^^^^ 
note the block sub-matrix of ^C-'") beginning with i*^ block 
row and column and ending with the j*^ block row and col- 
umn, while a[''J^^ will denote the [i,])^^'' block element of 
_4(i,m)^ Also X„ will denote a n X n block-diagonal matrix 
with diagonal elements /. 

The inverse of y^^^^^ is denoted by [y'^^-^^-^ = ^(^'^^ 
and satisfies the equation: 



y 



{1,N} ^{IM) 



I. 



N 



According to our notation we have = Q+ and G\ ^ 

G 



j^. The matrix J^f^'^) has the following structure: 

37(1^^) = ( ^"^^ ) , (B6) 



where = (0,0,...,-/) is a 1 x iV 
then write Eq. (B6I in the form 



1 block vector. We 



S+ 








. 


. 





. 


. 





. 





(Pl) 




u 



N 



G 



N.N 



In-1 
/ 



,(B7) 



where 



(B2) 



I \s a N' X N' unit matrix, and is a iV' x iV' matrix with 
all elements equal to zero. Similarly the matrices and V 
have the following block structure: 



M 



/Ml 
M2 

V 




V 




,G 



{1,N}T 
N-1,N 



] is a lxA^-1 block 



vector. From Eq (B7l we get the following four equations: 



y{l,N-l) gil-N) 



^N 



^{1,N-1) 



WnK=^n^i 
' CLN Uj] — , 



J^d'^-i) Um + Wat G^^;^) = , 



,(B3) 



I 



(B8) 



where denotes the diagonal mass-matrix for the ni = n 
layer and $ is a force-constant matrix whose off-diagonal 
terms correspond to coupling to sites within a layer. Hence 
the matrix = [—Muj'^ + V — S'l ~ 5^] has the following 
structure: 



Noting that [y^^-'^-'^'^]-'^ = g(i^^-i) we get, using the third 
equation above and the form of Wat: 



or G 



,N 



^i.N-1 ^N.N ' ror I 



l,2,...,iV- l.(B9) 



/ ai 
' -I 



[GY 



V 



0-2 

... 





-/ 
-/ 





a-N-i 



-/ 
a-N J 



where a/ = —Miuj"^ 



(B4) 



Now defin- 



ing Tl ji = Im\Ti^ ^ and with the form of S'^ ^ given in 



Eqs. (Bl l,(^B2|, we find that the expression for the transmis- 



(B5) 



equation. 



From the fourth equation in Eq. (B8I we get: 



G 



N-1,N 



"TV '-'N N 



I 



(BIO) 



sion coefficient reduces to the following form: 

TNiuj) = 4 Tr[TLiLo)G+{u)TR{L,)G],{Lo)] , 

where is the (IjN)*"^ block element of Q and G^ = 
[G^]^. We now show that G^ satisfies a simple recursion Setting 



We will now use Eqs. (|B9[ ),(B10l to obtain a recursion for 
G^^J^' = G% in Eq. \B5\], which is the main object of in- 
terest. Let us define P^-") = [Gi^n-i+ir^ where ^^C-™' = 
[•y(;,m)]-i^ Then setting i = 1 in Eq. (|B9]) and taking an 
inverse on both sides we get: 



N-l) 



(Bll) 



N 



1 in Eq. (B9l we get G 



1(1, ^v) 



tv-i'-^tv'tv^ using this in Eq. (|B10|) we get 
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[G 



N,N 



[gn - G 



Af-l,Ar-lJ 



Inserting this in the above 



equation we finally get our required recursion relation: 

The initial conditions for this recursion are: P(i'") = Im and 
p(i.i) = By proceeding similarly as before we can also 
obtain the following recursion relation: 



71+2, N) 



n=l,2,...,iV-]:^13) 



and P(i'^) can be recursively obtained using the initial con- 
ditions p(^+i.^) = Im and p(^'^) = ajy. Given the set 
{fli}, by iterating either of the above equations one can nu- 
merically find p(i'^) and then invert it to find G^^j^''. How- 
ever this scheme runs into accuracy problems since the nu- 
merical values of the matrix elements of the iterates grow 
rapidly. We describe now a different way of performing the 
recursion which turns out to be numerically more efficient. 
We first define 



r^, = p(i'^)[p(i'^-i)]-i 



From Eq. (B12l we immediately get: 

1 

rN = a-N 



rN-i 



(B14) 



(B15) 



system and reservoir contributions are separated. First we 
note that the form of the matrices ai for our specific problem 
is: ai = ci — iJ/^iSi — Si^n'^n where c/ = —Milu"^ + We de- 
fine system-dependent matrices by replacing 
ai,ajsi by Ci^c^ in the recursions for Ps'. Thus Q^^'"^ = 
p(i,»)(a^ ^ Ci,aN cn) and Q("'^' = p("'^)(ai ^ 
ci,ajv Cn). Clearly Qs' satisfy the same recursion as the 
Ps' with ai replaced by c/. Then using Eqs. (B12l,(B13l, 
and similar equations for the Qs' we get: 



p(l,JV) 

Qil,N) _ Qi2,N) _ Q(l,^-1) + Q(2,JV-1) 

Qil.N) _Q(2,N) \ / 1 

2(i,JV-i) _g(2,A'-i) M Si 



= (1 -Stv) 



(B17) 



From the recursion relations for the Qs' it is easy to see that 



^(1,^) _Q(2,W) 
Q(l.W-l) _Q(2,N-l) 

/ M Q(i^^-2) 



_q(2M-1) 
_0{2,N-2) 



,..Ti , 



where 



(B18) 



with the initial condition ri = ai. Then G^^'^^ is given by 



[P 



(B16) 



This form where at each stage ^ is evaluated turns out to 
be numerically more accurate. 

Finally we show that one can express G^^]^' in the form 
of a product of matrices. The product form is such that the 



T, = 



a-i -I 
I 



(B19) 



We then obtain Gj, = [P^^'^'>]-^. 

In our numerical calculations we use the recursion rela- 



tions in Eqs. (B15l,(B16) to evaluate the required Green's 



function. Computing the trace in Eq. (B5l then gives us the 



transmission coefficient as a function of frequency. 
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(a)Free boundaries (b)Fixed boundaries (c)Pinned lattice 



FIG. 1: A schematic diagram of a two-dimensional mass-disordered lattice of particles connected by harmonic springs and connected 
to heat baths at temperatures Tl and Tr. Red and green colours indicate particles of different masses. Pinning refers to the presence 
of a spring attaching a particle to the substrate. In (a) there is no pinning, in (b) boundary particles are pinned and in (c) all sites are 
pinned. 




FIG. 2: (color online) ID unpinned case with both free and fixed (INSET) boundary conditions: plot of the disorder averaged 
transmission T{uj) versus to for A = 0.4. The various curves (from top to bottom) correspond to lattices of sizes A'' = 64, 256, 1024 
respectively. 
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FIG. 3: (color online) Unpinned lattices with fixed BC in one direction and periodic in all others. 

Disorder averaged density of states obtained numerically from the eigenvalues of several disorder realizations in ID, 2D and 3D for 
lattice sizes A'' = 4096, 64, 16 respectively. Note that the low frequency behaviour is unaffected by disorder and one has oj*^"^ as a; — > 0. 
We set A = 0.8, k — 1 and averaged over 30 realizations in ID and over 10 realizations in 2D and 3D. In 2D and 3D there is not much 
variation in p{uj) for different disorder samples. Also shown are the density of states for the binary mass ordered lattices. 
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FIG. 4: (color online) Pinned lattices. 

Disorder averaged density of states obtained numerically from the eigenvalues of several disorder realizations in ID, 2D and 3D for 
lattice sizes A'^ = 4096, 64, 16 respectively. Note that low frequency modes are absent. We set k — 1, ko = 10.0 and A = 0.4 in 2D and 
A — 0.8 in ID, 3D. Averages were taken over 30 realizations in ID and 10 realizations in 2D, 3D. We find that in 2D and 3D there is 
not much variation in p{uj) for different disorder samples. Also shown are the density of states for the binary mass ordered lattices. 
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FIG. 5: (color online) 2D unpinned case with fixed BC for A = 0.95. (i) Plot of the disorder averaged transmission T(lo) versus lo. 
(ii) Plot of NT(lS). The range of frequencies for which T(lo) ~ IjN is indicated by the dashed line, (iii) Plot of p(tj) for binary mass 
ordered and single disordered sample, (iv) Plot of N^P~^ for single samples (smoothed data). We see that even though the allowed 
normal modes occur over a large frequency band « (0 — 12), transmission takes place in a small band « (0 — 1.25) and is negligible 



elsewhere. The IPR plots confirm that the non-transmitting states correspond to localized modes, 
decreasing with inrease of A^. 
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FIG. 6: (color online) 2D unpinned case with fixed BC for A = 0.8. 

TOP: Plot of the disorder averaged transmission T(tj) versus uj. The various curves (from top to bottom ) correspond to square lattices 
with A'' = 16,32,64,128,256,512 respectively. We see again that most modes are localized and transmission takes place over a small 
range of requencies. 

BOTTOM: Plot shows the IPR (P^^) as a function of normal mode-frequency tOp for the 2D lattice with A = 0.8. The curves are 
for TV = 16 (blue), 32 (green) and 64 (red). The inset plots N^P~^ and the collapse at low frequencies shows that these modes are 
extended. Also shown are two typical normal modes for one small (left) and one large value of for A'' = 64. 
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FIG. 7: (color online) 2D unpinned case with fixed BC for A = 0.2. 

TOP: Plot of the disorder averaged transmission versus ui. The upper-most curve corresponds to a binary-mass ordered lattice with 
A'^ = 128 while the remaining curves (from top to bottom) correspond to square lattices with A'^ = 16, 32, 64, 128, 256, 512 respectively. 
BOTTOM: Plot shows the IPR (P"^) and scaled IPR {N^P'^) as a function of normal mode- frequency uip. The curves are for A'^ = 16 
(blue), 32 (green) and 64 (red). 



17 
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FIG. 8: (color online) 2D unpinned lattice with fixed BC. 

Plot of disorder-averaged current J versus system size for different values of A. The error-bars show the actual standard deviations 
from sample-to-sample fluctuations. Numerical errors are much smaller. For A — 0.2, simulation data is also plotted. 
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FIG. 9: (color online) 2D unpinned case with fixed BC for A = 0.2. 

Plot of disorder-averaged temperature profile [Ti] for different system sizes obtained from simulations. 
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FIG. 10: (color online) 2D unpinned case with free BC for A = 0.8. 

TOP: Plot of the disorder averaged transmission T{uj) versus u. The various curves (from top to bottom) correspond to square lattices 
with TV = 16, 32, 64, 128, 256, 512 respectively. We see that transmission takes place in a small band ~ (0 — 2) of the full range ~ (0 — 6) 
of normal modes and as can be seen in the inset is negligible elsewhere. 

BOTTOM: Plot shows the IPR {P~^) as a function of normal mode-frequency ujp. The curves are for A'' = 16 (blue), 32 (green) and 
64 (red). In the inset we plot N^P~^ and the collapse at low frequencies shows that low frequency modes are extended. Also shown 
are two typical normal modes for one small (left) and one large value of P^^ for A'^ = 64. 




FIG. 11: (color online) 2D unpinned case with free BC for A 
Plot of the disorder averaged transmission Tiuo) versus ui for, 
respectively. Note the linear form at smah uj. 



= 0.2. 

■. The curves (from top to bottom) are 



for N = 16,32,64,128,256,512 
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FIG. 12: (color online) 2D unpinned case with free BC. 

Plot of disorder-averaged current J versus system size for two different values of A. The error- bars show standard devations due to 
sample-to-sample fluctuations. Numerical errors are much smaller. 




CO 



FIG. 13: Plot of the effective mean-free path lc« — NT{uj)/uj''~^ in 2D with A = 0.8. The insets show £cS for the ordered lattices 
with a single mass. An u)~^ behaviour is observed in a small part of the diffusive region. The fixed BC data is highly oscillatory and 
has been smoothed. 
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FIG. 14: (color online) 2D pinned case for A — 0.4 and ko — 10.0. 

TOP: Plot of the disorder averaged transmission T{uj) versus uj. The various curves (from top to bottom) correspond to lattices with 
A'' = 16, 32, 64 respectively. Here we choose 7 — \/IO. 

BOTTOM: Plot of the IPR (P~^) as a function of normal mode-frequency ujp. The curves are for = 16 (blue), 32 (green) and 64 
(red). Also shown are two typical normal modes for one small (left) and one large value of for TV = 64. 
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FIG. 15: (color online) 2D pinned case for A = 0.4 and ko — 2.0. 

Plot of the disorder averaged transmission T{uj) versus ui . The various curves (from top to bottom) are for A'^ = 16, 32, 64, 128, 256, 512 
respectively. Here we choose 7 — \/2. 




FIG. 16: (color online) 2D pinned case for A — 0.4. 

Plot of disorder-averaged current J versus system size for two different values of ko- Error bars show standard deviation due to disorder 
and numerical errors are much smaller. Note that the standard deviation do not decrease with system size for higher ko. 
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FIG. 17: (color online) 2D pinned case for A = 0.4 and ko = 10.0. 

Plot of disorder-averaged temperature profile [Ti] for different system sizes. The plots are from simulations and here we choose 7 = vTO. 
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FIG. 18: (color online) 3D unpinned case with fixed BC for A — 0.8. 

TOP: Plot of the disorder averaged transmission T{lu) versus uj. The inset shows the same data multiplied by a factor of A'^. 
BOTTOM: Plot of the IPR (-P~^) and scaled IPR {N''^P~^) as a function of normal mode-frequency LUp for a fixed disorder-realization. 
The curves are for N — 8 (green) and 16 (red). 



25 



3^ 



— N = 


8 


-- N = 


16 


.... N = 


32 




FIG. 19: (color online) 3D unpinned case with fixed BC for A — 0.2. 

TOP: Plot of the disorder averaged transmission T{uj) versus uj. The uppermost curve is the transmission curve for the binary mass 
ordered lattice for A'' = 16. 

BOTTOM: Plot of IPR (P^^) and scaled IPR (N'^P^^) as a function of normal mode-frequency LUp for a fixed disorder-reahzation. 
The curves are for A'^ = 8 (green) and 16 (red). 




FIG. 20: (color online) 3-D unpinned case with fixed and free BCs. 

Plot of disorder-averaged current J versus system size for two different values of A. The data for A — 0.2 is from simulations. The 
error-bars show standard deviations due to disorder and numerical errors are smaller. 
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FIG. 21: (color online) 3-D unpinned case with fixed BC for A = 0.2. 

Plot of disorder-averaged current density J (with the definition J — I/N2) versus N2/N for different fixed values of A''. We see that 
the 3D limiting value is reached at quite small values of /N. 
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FIG. 22: (color online) 3-D unpinned case with fixed BC for A = 0.2. 

Plot of temperature profile Ti in a single disorder realization for different system sizes. The plots are from simulations. 
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FIG. 24: Plot of the effective mean-free path les = NT{uj) /uj'''~^ in 3D with A = 0.8 for fixed and free BCs. The insets show 4fr 
for the ordered system with a single mass. An a;~* behaviour is observed in a small part of the diffusive region. The fixed BC data is 
highly oscillatory and has been smoothed. 
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FIG. 25: (color online) 3-D pinned case for A — 0.2 and ko — 10.0. 
TOP: Plot of the disorder averaged transmission T(uj) versus uj. 

BOTTOM: Plot of the IPR (P"^) and scaled IPR (TV^P"^) as a function of normal mode-frequency ujp. The curves are for iV 
(green) and 16 (red). 
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FIG. 26: (color online) 3-D pinned case for A — 0.8 and ko — 10.0. 
TOP: Plot of the disorder averaged transmission T{lo) versus ui. 

BOTTOM: Plot of the IPR (P"^) scaled IPR {N^P''^) as a function of normal mode-frequency ujp. The curves are for N 
and 16 (red). 
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FIG. 27: (color online) 3-D pinned case. 

Plot of disorder-averaged current J versus system size for different values of ko and A. The data sets for A 
of ko are from simulations while the data for A = 0.8 is from numerics. 
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FIG. 28: (color online) 3D pinned case for A = 0.2 and fco = 10.0. 

Plot of temperature profile Ti in a single disorder realization for difi'erent system sizes. 



The plots are 



from simulations. 



